%% Reflectivity Analysis
% Based on Effective Density Model
% May-26-2016
% use plcol.m
close all
clear all
clc

%FileName='sph01.dat';
FileName='sph02.dat';
TitleName='Sphingomylin';
ScanNumber=[
    112,113;   %  #1 sph01.dat SPM + 100 mM KCl; pH 3;
    140,141;   %  #2 sph01.dat SPM + 100 mM KCl; pH 4;
    177,178;   %  #3 sph01.dat SPM/Fe + 100 mM KCl; pH 3;
    251,252;   %  $4 sph01.dat SPM/Fe +100 mM KCl; pH 7.08
    321,322;   %  $5 *  sph01.dat SPM/Fe +100 mM KCl; pH 4
    27,28;     %  $6 sph02.dat SPM/Fe +100 mM KCl; pH 2
    92,93;     %  #7 sph02.dat SPM/Fe +100 mM KCl; pH 5
    151,152;   %  #8 sph02.dat SPM/Fe +100 mM KCl; pH 5
    215,216;   %  #9 sph02.dat SPM/Fe +100 mM KCl; pH 3
    274,275;   %  #10 sph02.dat SPM/Fe +100 mM KCl; pH neutral
    333,334;   %  #11 sph02.dat SPM/Fe +100 mM KCl; pH neutral
    341,342;   %  #12 sph02.dat SPM/Fe +100 mM KCl; pH neutral
    ]
RefTxtFile={
    'sph01_112_ref.txt';   %  #1 sph01.dat SPM+ 100 mM KCl; pH 3;
    'sph01_141_ref.txt';   %  #2 sph01.dat SPM+ 100 mM KCl; pH 4;
    'sph01_177_ref.txt';   %  #3 sph01.dat SPM/Fe+ 100 mM KCl; pH 3;
    'sph01_251_ref.txt';   %  #4 sph01.dat SPM/Fe+ 100 mM KCl; pH 7.08;
    'sph01_321_ref.txt';   %  #5 sph01.dat SPM/Fe+ 100 mM KCl; pH 4;
    'sph02_27_ref.txt';   %   #6 sph02.dat SPM/Fe+ 100 mM KCl; pH 2;
    'sph02_92_ref.txt';   %   #7 sph02.dat SPM/Fe+ 100 mM KCl; pH 5;
    'sph02_152_ref.txt';   %   #8 sph02.dat SPM/Fe+ 100 mM KCl; pH 5;
    'sph02_215_ref.txt';   %   #9 sph02.dat SPM/Fe+ 100 mM KCl; pH 3;
    'sph02_274_ref.txt';   %   #10 sph02.dat SPM/Fe+ 100 mM KCl; pH neutral 
    'sph02_333_ref.txt';   %   #11 sph02.dat SPM/Fe+ 100 mM KCl; pH neutral 
    'sph02_341_ref.txt';   %   #12 sph02.dat SPM/Fe+ 100 mM KCl; pH neutral 
    
    }
ParamToSave={
    'sph01_112_ref_par.mat';   %  #1  sph01.dat SPM+ 100 mM KCl; pH 3
    'sph01_141_ref_par.mat';   %  #2  sph01.dat SPM+ 100 mM KCl; pH 4
    'sph01_177_ref_par.mat';   %  #3  sph01.dat SPM/Fe + 100 mM KCl; pH 3
    'sph01_251_ref_par.mat';   %  #4  sph01.dat SPM/Fe + 100 mM KCl; pH 7.08
    'sph01_321_ref_par.mat';   %  #5  sph01.dat SPM/Fe + 100 mM KCl; pH 4  
    'sph02_27_ref_par.mat';   %   #6  sph02.dat SPM/Fe + 100 mM KCl; pH 2
    'sph02_92_ref_par.mat';   %   #7  sph02.dat SPM/Fe + 100 mM KCl; pH 5
    'sph02_152_ref_par.mat';   %   #8  sph02.dat SPM/Fe + 100 mM KCl; pH 5
    'sph02_215_ref_par.mat';   %   #9  sph02.dat SPM/Fe + 100 mM KCl; pH 3
    'sph02_274_ref_par.mat';   %   #10  sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    'sph02_333_ref_par.mat';   %   #11  sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    'sph02_341_ref_par.mat';   %   #12  sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    }

OptimParamSave={
    'sph01_112_ref_optim.txt';   %  #1  sph01.dat SPM+ 100 mM KCl; pH 3
    'sph01_141_ref_optim.txt';   %  #2  sph01.dat SPM+ 100 mM KCl; pH 4
    'sph01_177_ref_optim.txt';   %  #3  sph01.dat SPM/Fe + 100 mM KCl; pH 3
    'sph01_251_ref_optim.txt';   %  #4  sph01.dat SPM/Fe + 100 mM KCl; pH 7.08
    'sph01_321_ref_optim.txt';   %  #5  sph01.dat SPM/Fe + 100 mM KCl; pH 4
    'sph02_27_ref_optim.txt';    %  #6  sph02.dat SPM/Fe + 100 mM KCl; pH 2
    'sph02_92_ref_optim.txt';    %  #7  sph02.dat SPM/Fe + 100 mM KCl; pH 2
    'sph02_152_ref_optim.txt';   %   #8  sph02.dat SPM/Fe + 100 mM KCl; pH 5
    'sph02_215_ref_optim.txt';   %   #9  sph02.dat SPM/Fe + 100 mM KCl; pH 3
    'sph02_274_ref_optim.txt';   %   #10 *  sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    'sph02_333_ref_optim.txt';   %   #11 * sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    'sph02_341_ref_optim.txt';   %   #12  sph02.dat SPM/Fe + 100 mM KCl; pH neutral
    }

RefDataSelect=12;
data1=plcol(FileName,ScanNumber(RefDataSelect,1),'q1','ref','sref');
data2=plcol(FileName,ScanNumber(RefDataSelect,2),'q1','ref','sref');
q=[data1(:,1);data2(:,1)];
ref=[data1(:,2);data2(:,2)];
ref_max=max(ref);
sref=[data1(:,3);data2(:,3)];
H=ploterr(q,ref,[],sref,'ko','hhy',0.1);
set(H,'LineWidth',3,'MarkerSize',10,'MarkerFaceColor','w');
set(gca,'yscale','log');
xdata=q;
ydata=ref;
ydata_err=sref;
Txt_Save(RefTxtFile{RefDataSelect},xdata,ydata,ydata_err);

%% %%--------------------------------
P_length=16;
C_length=10;
C=zeros(1,C_length);
C(1)=100; % 100 slices
C(6)=0;
C(5)=6.8285; % 2*pi/lambda lambda=0.9201 Angstrom
%% %%--------------------------
smlnum=1e-6;
clc
offset_min=-0.002;offset_max=+0.002;               % P[1]

P2_min=-smlnum;P2_max=smlnum;                          % P[2]

I0_min=0.8*ref_max; I0_max=1.2*ref_max;                           % P[3]
ED_sub_min=0.334-smlnum;ED_sub_max=0.334+smlnum;      % P[4]
mu_sub_min=0;mu_sub_max=smlnum;                     % P[5]

d1_min=3; d1_max=30;                             % P[6]
rho1_min=0.334;rho1_max=1;              % P[7]
sigma1_min=3; sigma1_max=20;              % P[8]

d2_min=3; d2_max=20;                             % P[9]
rho2_min=0.317;rho2_max=1;                      % P[10]
sigma2_min=1; sigma2_max=10;                      % P[11]

d3_min=3; d3_max=20;                             % P[12]
rho3_min=0.317;rho3_max=1;                      % P[13]
sigma3_min=1; sigma3_max=10;                      % P[14]

sigma4_min=1; sigma4_max=10;                     % P[15]

P16_min=-smlnum;P16_max=smlnum;                % P[16]

LB=[
    offset_min;P2_min;I0_min;ED_sub_min;mu_sub_min;...
    d1_min;rho1_min;sigma1_min;...
    d2_min;rho2_min;sigma2_min;...
    d3_min;rho3_min;sigma3_min;...
    sigma4_min;...
    P16_min
    ];

UB=[
    offset_max;P2_max;I0_max;ED_sub_max;mu_sub_max;...
    d1_max;rho1_max;sigma1_max;...
    d2_max;rho2_max;sigma2_max;...
    d3_max;rho3_max;sigma3_max;...
    sigma4_max;...
    P16_max
    ];

%%
fval_min=1e7;
M=200; % Number of shakeup
coef=rand(7,M);% 4 is the number of parameter to shake up
Param_Opt=0.5*(LB+UB);
%Param_Opt(3)=1;
Param_stack=zeros(17,M); % Last one is chisq
for k=1:M
    ptry=Param_Opt;
    ptry(9:15)=coef(:,k).*(UB(9:15)-LB(9:15))+LB(9:15);
    [x1,fval,residual,exitflag,output]=optim_fitref_ED_3box(ptry,C,xdata(xdata<0.6),ydata(xdata<0.6),ydata_err(xdata<0.6),LB,UB);
    if fval<fval_min
        fval_min=fval;
        Param_Optim=x1;
    end
   Param_stack(1:16,k)=x1;
    %chisq_reduce=fval/(length(xdata));
    
   Param_stack(17,k)=fval;
    disp(strcat('loop :',num2str(k),' chisq = ',num2str(fval)))
end
save(ParamToSave{RefDataSelect},'Param_stack');
Txt_Save(OptimParamSave{RefDataSelect},Param_Optim);
%%
close all
xfit=[min(xdata):1e-4:max(xdata)]';
figure

hold on
C(2)=0;
qcutoff=Param_Optim(16);
dq=Param_Optim(1);

rfresnel=FresnelReflection(xdata-dq,0.0217,0,1.54);
ydata_norm=ydata./(rfresnel);
H=ploterr(xdata,ydata_norm,[],ydata_err./rfresnel,'ko','hhy',0.1)
set(H,'LineWidth',3,'MarkerFaceColor','w','MarkerSize',12);
rfresnel_fit=FresnelReflection(xfit-dq,0.0217,0,1.54);
yfit=fitref_ED_3box(Param_Optim,C,xfit)./rfresnel_fit;
plot(xfit,yfit,'b-','LineWidth',3);
hold off
set(gca,'yscale','log');
set(gca,'xlim',[0,0.55])
grid on
axis square

